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ABSTRACT 

We present a new approach to gravitational lens massmap reconstruction. Our massmap solutions 
perfectly reproduce the positions, fluxes, and shears of all multiple images. And each massmap 
accurately recovers the underlying mass distribution to a resolution limited by the number of multiple 
images detected. We demonstrate our technique given a mock galaxy cluster similar to Abell 1689 
which gravitationally lenses 19 mock background galaxies to produce 93 multiple images. We also 
explore cases in which as few as four multiple images are observed. Massmap solutions are never 
unique, and our method makes it possible to explore an extremely flexible range of physical (and 
unphysical) solutions, all of which perfectly reproduce the data given. Each reconfiguration of the 
source galaxies produces a new massmap solution. An optimization routine is provided to find those 
source positions (and redshifts, within uncertainties) which produce the "most physical" massmap 
solution, according to a new figure of merit developed here. Our method imposes no assumptions 
about the slope of the radial profile nor mass following light. But unlike "non-parametric" grid-based 
methods, the number of free parameters we solve for is only as many as the number of observable 
constraints (or slightly greater if fluxes are constrained). For each set of source positions and redshifts, 
massmap solutions are obtained "instantly" via direct matrix inversion by smoothly interpolating 
the deflection field using a recently developed mathematical technique. Our LensPerfect software is 
straightforward and easy to use and is made publicly available via our website. 

Subject headings: gravitational lensing — methods: data analysis — galaxies: clusters: general — 
cosmology: dark matter 
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1. INTRODUCTION 

Simulations of structure formation in a ACDM uni- 
verse have provided concrete predictions for the form 
of Da rk Matter halos over a wid e rang e of scales (e.^., 
Merritt et al.^ l2QQ6l : iBett et all 120071 : iDiemand etall 
20071 ). These mass distributions, quantified in terms of 
radial profile, ellipticity, and level of substructure, are 
among the key predictions of ACDM theory. Uncertain- 
ties do persist, especially with regard to baryons, absent 
from most Dark Matter simulations. Baryons are found 
to alter the inner profiles and ellipticities of halos, es- 
pecially on galaxy scales (e.g., IKa zantzid is et al.l l2004l : 
iGustafsson et all l2006l : button et al. 2007). 

Testing these predictions in detail observationally has 
proven challenging. Gravitational lensing provides us 
with our most direct tool for mapping the distributions of 
mass (predominantly Dark Matter) within and surround- 
ing galaxies and galaxy clusters. But massmaps recov- 
ered by this method are of much lower resolution than 
simulated Dark Matter halos. Improvements in imaging 
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quality both from the ground and space now allow for a 
more definitive measurement of lensing effects, motivat- 
ing new techniques to take full advantage of this advance. 

Deep multi-wavelength high-resolution images of five 
galaxy clusters (Abeh 1689, Abeh 2218, Abeh 1703, 
MS1358, and CL0024) have been obtained with A CS, 
the Advanced Camera for Surveys (Ford et aLll2003l) on- 
board the Hubble Space Telescope. The deepest of these, 
taken of A1689, has revealed an unprecedented number 
of multiple images, with 106 images of 30 s ource galax- 
ies identified in the original analysis by Broadhurst e t al.l 
(|2Q05) . This should allow one to reconstruct a relatively 
high resolution massmap of the cluster's Dark Matter 
halo. 

iBroadhurst et al.l (|2005l ) used a novel massmap param- 
eterization to reproduce the 106 observed multiple im- 
ages well but not exactly (Fig. [1^), with positional off- 
sets of 3^2 RMS on average in the image plane. Note 
that these are very significant offsets, as image posi- 
tions can typically be determined to 0V05 (the resolu- 
tion of ACS pixels) or better. Subsequent analyses us- 
ing similar methods (with different parameterizations for 
the massmap and revised multiple image lists) improved 
only mar ginally on these po s itional offsets to 2^5 RMS 
at best dZekser et al.l 120061 : iHalkola et all l2006l . l2007l : 
iLimousin et al.ll2006 ^. It is possible that cosmic variation 
of mass density along the multiple sight lines is respon- 
sible for some of this scatter. But it is more likely that 
line of sight variations are an additional source of error 
(present in any lensing analysis) on top of the models' 
inability to properly reproduce all the image positions. 

The aforementioned analyses all employ extensions 
of techniques developed for use with far fewer multi- 
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Fig. 1. — Left: Multiple images of object #1 from (Broadhurst et al.''2005") delensed back to the source plane using their best fit defiection 
field. The scatter in the source plane (~ 2") is typical of all model-based massmap reconstruction methods used to date. Right: The same 
multiple image positions all delensed back to a single source position. These defiections can be fit exactly by a LensPerfect solution. 



pie images. A simple parametric model consisting of 
an elliptical mass profile plus external shear can ac- 
curately reprodu ce the positions of four multiple im- 
ages (e.g., Keeton et al.l fl997). But such models strug- 
gle when confronted with additional multiple images. 
Multi-component ("halo" -h "galaxy") models are re- 
quired to obtain d e cent fits to galaxy cluste r lensing (e.g., 
iKneib et al.lll996l : iBroadhurst et al.ll2QQQ[ ). The A1689 
multi-component models discussed above don't contain 
enough free parameters to perfectly fit the data. Yet 
they already contain so many parameters that navigat- 
ing the parameter space proves challenging. It may be 
impossible to find the best solution with so many free 
parameters given our current computational capabilities. 
Advances in computing power may someday allow para- 
metric models the freedom necessary to produce perfect 
fits. Allowing galaxy subhalos to drift in position and 
vary individually in mass could dramatically improve the 
fits while breaking free from the assumption that Dark 
Matter subhalos strictly follow the light distribution. 

The degree to which Dark Matter follows light (and/or 
gas) is an important outstanding question. The exciting 
discoveries in this area are currently being made by weak 
gravitational lensing analyses which make no assump- 
tions about the underlying mass distributions. The "Bul- 
let Cluster" finding (Clowe et al. 2006) that gas has been 
stripped from two colliding Dark Matter halos would not 
have been possible had the authors assumed mass follows 
light from the outset. A similar cluster collision observed 
along the line of sight appears to have left a Dark Mat- 
ter ring around CL002 4. This ring wa s also detected by 
weak lensing analysis (|Jee et al.ll2007[ ). 

For many years now, we have obtained assumption- 
free massmap reconstructions fron i direct inversion 
of weak gravitational lensing data (jTyson et al.l 1 19901 : 
iKaiser fc Squires 1993). Here we present the first 
method to do the same given strong gravitational lens- 
ing data (multiple images). Our method is not entirely 
assumption free, as a few basic assumptions about the 
distribution of mass can be helpful in selecting the most 
"physical" among the possible massmap solutions (§ 12. 4p . 

Model-free massmap reconstructions have previously 
been obtained for Abell 1689 using strong lensing data 
(though not via direct inversion) . As these massmaps are 
more fiexible than model-based solutions, they should 



produce better fits to the data. But this promise has 
yet to be fully realized. The SLAP method (|Diego et al.l 
'2005a) computes massmaps on an adaptive mesh grid, 
and fits the data to a desired level of scatter. But 
when this level is set too low, the solutions become "bi- 
ased" with "a lot of substructure". Their best solution 
for A 1689 leaves sc atters of ~ 3^' in the source plane 
(jDiego et al.l l2005bl: Diego 2007, private communica- 
tion). Meanwhile . Is aha et al.l (|2006l ) use a method called 
PixeLens to obtain pixel-based (fixed grid) massmaps 
which perfectly reproduce all multiple image positions. 
But computational constraints prevent them from fitting 
more than 30 images at a time. The massmap we obtain 
for A1689 (to be presented in an upcoming paper) per- 
fectly reproduces the positions of 135 multiple images (as 
in Fig.[T]3) and thus has about twice the resolution as the 
PixeLens massmap in each spatial direction (as dictated 
by the density of multiple images). 

Massmaps may be further improved by incorporating 
constraints beyond simple image positions. Images which 
are resolved and exten ded should be properly reproduced 
by the mass jnodejjjmrren fc Dvd l2003l : [S uvu et al.l 
120061 : iKoopmans et al.[r2006[ ). 

Even unresolved images yield the information encoded 
in their fiuxes. A simple mass model which accurately 
reproduces "quad" image positions may or may not ac- 
curately reproduce the image fiux ratios. Discrepant 
cases are known as "fiux anomalies" . If other causes (mi- 
crolensing, time delays, and dust extinction) can be ruled 
out or accounted for, then small substructure (~ lO^M© 
subhalos) is generally invoked as the most likely expla- 
nation for an observed anomaly. 

But this explanation has perhaps been invoked too of- 
ten. The amount of substructure observed in simulations 
may not be sufficient t o produce fiux anomalies as of- 
ten as observed (M etc alf et al.l[2003: lAmara et all l2006l : 
iMaccio Mirandal l2006l : iDiemand etld] l2007f ). One 
way to resolve this possible discrepancy is to obtain 
macro mass model solutions which reproduce the ob- 
served fiuxes without resorting to smaller substructure. 

To address this issue, Evans & Witt (2003) developed 
a direct inversion massmap reconstruction method ca- 
pable of perfectly fitting the observed positions and 
fiuxes of four multiple images. While providing rea- 
sonable solutions for the lensed systems Q2237+0305 
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and PG1115+080, their solution for B1422+231 was 
clearly unphysical. While it is unclear how exactly to 
quantify a massmap's physicality, the authors charac- 
terized their B1422+213 model as simply too "wiggly" 
to be plausible. D eveloping the method a bit further, 
ICongdon fc Keeton l (12005) produced a more reasonable 
solution for B1422+231, but were less successful with 
B2045+265 and B1933+503. For the unsuccessful cases, 
the authors argue that small scale substructure is the 
most likely explanation for the observed flux anomalies. 
But perhaps their models were simply not flexible enough 
to obtain physical solutions for these systems. We will 
revisit this question in an upcoming paper. 

Our new LensPerfect method uses direct inversion to 
obtain assumption-free massmap solutions which per- 
fectly reproduce all multiple image positions. Multiple 
image knots and fluxes may also be perfectly constrained. 
LensPerfect is made possible by a recent advance in the 
field of mathematics. The essential tool is a method that 
produces a curl-free interpolation of vectors given at scat- 
tered data points. A set of observed multiple image posi- 
tions and redshifts along with assumed source positions 
defines a deflection field at the image positions. Once we 
interpolate this deflection field across the entire field, we 
take the gradient, multiply by 2, and "instantly" obtain 
our perfect massmap solution. 

The interpolation of data given at scattered points 
is a complex problem without a unique solution. One 
method of attacking such problems involves the use of 
radial basis functions (RBFs). RBFs have been used to 
interpolate scattered data since the early 1970's (Hardv 
EMI) and are used in many applications today. By 
the 1980's, RBFs had been applied to interpolate vec- 
tor fields. And in 1994, a method was developed capa- 
ble of yielding divergence-free interpolations of scattered 
vectors (Narcowich & Ward 1994). Finally, the curl- free 
analo g of this method was developed last year (Fuselier 
I2006L 120071 Here we apply this new method to gravita- 
tional lensing analysis. 

We describe our method along with the necessary grav- 
itational lensing equations in ^ In § [31 we demonstrate 
applications of the method, including the recovery of a 
known massmap given 93 multiple images (§ 13. ip . We 
also explore the solutions obtained when only a hand- 
ful of multiple images are available (§ 13. 4p . And we 
demonstrate the gains made by constraining extra im- 
age knots in § 13.61 The method for adding flux and/or 
shear constraints is given in § 12.61 In § |4] we discuss the 
relative merits of model-based and model-free methods, 
along with the potential of a hybrid method among other 
techniques that may become possible with LensPerfect. 
Finally, we provide a summary in § [5l The LensPerfect 
software and more information are available at our web- 
site.^ 

2. METHOD 

2.1. Image Deflection 

Image deflection by a gravitatio nal lens is go verned 
by a few simple equations (e.g., IWambsganssI [1998). 
The relativistic bending of light due to a mass M at a 
distance R away is twice that expected from Newtonian 
physics: a = AGM/c^R given Newton's gravitational 

^ http://www.its.caltech.edu/%7Ecoe/LensPerfect/ 




Fig. 2. — Light ray deflection by a gravitational lens of mass 
n. In the absence of k,^ the galaxy would appear at its true, or 
"source", position (3. The intervening mass deflects its light by an 
amount a to position 6. The deflection angle a on our sky is related 
to the actual bend angle a. of the light ray via olDs = olDls- Ds-, 
Dls^ and Dl are measured as angular diameter distances. 

constant G, and the speed of light c. In a gravitational 
lens, it is generally safe to assume that all of the 
deflection occurs in the plane of the lens. (This is known 
as the thin lens approximation.) Given the projected 

mass surface density distribution i<i{0) of the lens, we 
can derive the deflection of light d{0) = — [3 from its 
true position on the sky jS to that at which we observe 
it e (Fig. ED: 



2 ' 



(1) 



with the corresponding less-intimidating inverse relation: 



(2) 



The surface density n = T^/T^crit is defined in units of 
the critical density at the epoch of the lens. The critical 
density is that generally required for multiple images to 
be produced. It is a function of source redshift as given 
by: 



-'crit 



Ds 



4nG DlDls' 



(3) 



involving a ratio of the angular-diameter distances from 
observer to source Ds = Da{0^zs)^ observer to lens 
Dl = Da{0^zl), and lens to source Dls = Da{zl, zs). 
For a flat universe = l^rn + ^A = 1), angular- diameter 
distances are calculated as follows (Fukugit a et al.lfl992L 
filled beam approximation; see also Hogg 1999h: 



Da{zi,Z2) 



dz' 



1 + Z24 H{z'y 



(4) 
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where the Hubble parameter varies with redshift as: 



(5) 



Thus the critical density T^crit is a function of the 
source redshift. This follows because the deflection angle 
a is a function of source redshift. As source redshift 
decreases, the light bend angle a remains constant, 
which (imagine moving the galaxy in Fig. [2] inward 
along the top blue arrow) requires the image deflection 
to decrease by the distance ratio: 



(6) 



where doo is the deflection for a source at infinite redshift. 

Thus the problem of massmap reconstruction can 
be reduced to determining the deflection field with all 
deflections scaled to a common redshift (e.g., doo)^ at 
which point we simply take the divergence and divide 
by 2 to obtain the massmap (Eq. [2]). The deflection 

field d{6) = 6 — P may be measured at the multiple 

image positions once source positions /3 are determined 
(see § 12. 4p . However, in order to take its divergence, 
the deflection field must be solved for as a continuous 
function of position (or at least defined on a regular 
grid). Our interpolated deflection field must also be 
curl-free: 



V X a = 0. 



(7) 



This follows from Equation [T] in the same way that 
we find an electric field due to a static distribution of 
charge is also curl-free. One way to demonstrate this 
is to use the substitution Vln^ = 0/0'^ to define the 
lensing potential: 



such that 



iIj{0) = - I (fe'i^{e') In 

TT 



d — vv^. 



(8) 



(9) 



The fact that the deflection field d{0) may be written as 

the gradient of a scalar field il^iO) guarantees that it has 
no curl. 

2.2. Curl-Free Vector Interpolation 

iFuselierl (|2QQ6l , I2QQ7 ) has shown how to obtain a curl- 
free interpolation of a vector field given on scattered data 
points. The method is general enough to be applied to 
an arbitrary number of dimensions, but will be discussed 
here for the 2-D case. The interpolated vector field is 
constructed using Radial Basis Functions, or RBFs. An 
RBF is a positive definite function with radial symmetry 
(j){R/Ro). The fact that it is positive definite^ guarantees 

^ An m X m matrix- valued function (f) is positive definite on M"^ 
if given any finite, distinct set of points X := {xi, . . . ,xn} C M"' 

— Xk)oik > for all q;i, Qiiv in M"^. 

j,k 



we have aj 



that our interpolation matrix (b elow) will have a solution 
(jWendlandl 120051 : lFuselied[2006h . The scale length Ro is 
input by the user, and as we show in §3.4[ its freedom is 
closely related to the classic "mass-sheet" degeneracy. 

In our case we must choose smooth RBFs that are at 
least three times continuously different iable, or "C^", 
this to ensure our final massmap has no discontinuities. 
Then a set of two curl-free basis vectors are constructed 
at each data point (image position) by simply taking 
derivatives of the basis function: 



Vb = (p,yxX + (p,yyy, 



(10) 

(11) 



where (j)^xx is the 2nd derivative of (j){R/Ro) with respect 
to X, etc. The vector at that data point is given as a 
sum over all data points (with coefficients to be solved 
for): 



i 

Rewriting this in matrix form we have: 



(12) 



For example. 
d{02): 



H = [<P,i.v)] [c] ■ (13) 
in the case of 2 data points, d{Oi) and 
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\xy21 4^,yy2i 



^,XXi2 Y,yxi2 
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^a2 




_ ^b2 _ 



where, e.g.: 



dydx 



(14) 



(15) 



R=Ri2 



with the derivative being evaluated at = R12 = 
. We may now simply solve for the coefficients 
via the matrix inverse: 



M = [«] • (16) 

We may guarantee that this matrix inversion yields 
a solution by selecting a positive definite function for 
(j). Wendland functions are especially suitable choices, 
and here we use this "C^", or 6 times continuously 
different iable, Wendland function known as 1^7,3 (^): 



(1 
0, 



|C|)«(32|C|3 + 25|C|2 + 8|C| 



1), ICI<1 
ICI>1 
(17) 



where ^ = R/Rq. The function is very similar to 
a Gaussian with a peak of 1 and FWHM ~ 0.50i?o 

(cr - 0.21i?o) 
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In principle, a Gaussian could be used in place of ^7,3. How- 
ever it is well known (to mathematicians) that the interpolation 
matrix is ill-conditioned when a Gaussian is used. Further, it is 
easier for a computer to evaluate a polynomial than a Gaussian. 
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Fig. 3. — Massmap basis functions constructed as derivatives 
of the Wendland function Wr^s. a) ^Ca, b) ^c^, c) {na + Kb)/V^- 
d) Profile cut through the "major axis". Note that outside the 
domain [-1, 1], the function equals zero. 

Of relevance in this work are derivatives of this func- 
tion which serve as basis functions for our massmaps. 
Having solved for the deflection field [a] = [(l)^(xy)] H , 
the massmap may then be calculated analytically as: 



N 



(18) 



— 2 ^ V [^ai{4^,xxxi H~ 4^,xyyi) H~ (^bi{4^,yxxi H~ 4^,yyyi)^^^) 

i 

The basis functions: 



2 \Y.,xxx 



',yxx 



jxyyj 

\yyy) 



(20) 
(21) 



are shown in Figure [3l Different combinations of these 
two basis functions simply serve to rotate it about its 
axis and change its amplitude. 

The basis functions, with amplitudes and orientations 
solved for, are placed at the positions of the multiple im- 
ages. The resulting coefficients are generally many orders 
of magnitude greater than the amplitude of the massmap. 
These large mass components all cancel out nearly per- 
fectly with the small "residuals" being the massmap so- 
lution. Together, these mass components form a very 
flexible basis set. Although the basis functions are de- 
rived from radial basis functions, no symmetry conditions 
(radial, or otherwise) are imposed on the total massmap 
solution. 

We note that other RBFs may be used in place of 
VF7,3, including other Wendland functions, thin plate 
splines, multiquadrics, power laws, and even Gaussians 
(e.g., Fuselieri l2QQ6[ ). Wendland functions are especially 
well-suited for use in the interpolation matrix, but other 
functions may be explored in future work. 

2.3. Other Observables 

For reference, we provide here expressions for other 
observables as functions of our basis function. The lens- 
ing potential (Eq. [8]) is defined as a sum of derivatives 
of our basis function: 



Other observables may then be derived as derivatives of 
this potential. The deflection field is given as the gradient 
of the potential: 



We also obtain surface mass density 

^ ~ \ i'^jxx + ^^.yy) 1 



shear 



7+ = H^,^^ -^,2/2/)' 
Tx — ,xyi 



2 2,2 
7 =7++7x: 



the inverse of the magnification 



= 1 - i^.xx 



■7 

,2 



and time delays 



Ar = 



(1 + ^l)Ds 
cDlDls 



(23) 
(24) 



(25) 
(26) 
(27) 



(28) 
(29) 



(30) 



^ = Ca(t>,x + Cb4>,y 



(22) 



2.4. Source Positions and Massmap Physicality 

Image deflections are defined by the combination of im- 
age positions, source positions, and redshifts. But source 
positions are not known in practice and redshift measure- 
ments may contain large uncertainties. In our method, 
these become our free parameters. Each set of source 
positions and redshifts will redefine the deflection field, 
yielding a different LensPerfect massmap solution. We 
iterate to find the most plausible solutions including a 
single "best" solution, as we describe below. An ideal 
optimization procedure would lead us directly to the true 
source positions and redshifts. The optimization proce- 
dure we have developed may not be ideal, but we believe 
it produces an accurate "best" solution (see § 13. ip . sim- 
ilar to that obtained with the true source positions. We 
describe our optimization procedure in the next subsec- 
tion (§ 12. 5p after first describing our method for rating 
solutions. 

LensPerfect will return a perfect massmap solution for 
any set of input source positions and redshifts. But only 
certain sets will yield physical solutions: those with pos- 
itive mass everywhere within the convex hull, where our 
solutions are constrained. And some sets will yield "more 
physical" solutions than others. An extreme example of a 
"less physical" solution would be a "donut" solution com- 
prised of a high density ring surrounding a lower density 
center. By developing a method to rate different solu- 
tions as more or less physical, we can iterate over possible 
source positions to find those which yield plausible solu- 
tions, including the most plausible or "most physical" 
massmap solution. Our goal is to discard unreasonable 
solutions without biasing our result toward our concept 
of what a massmap should look like (that it should have 
an NEW profile, for example). 
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In the PixeLens method, a series of rigid criteria 
are used t o disting:uish ^^physica l" from "unphysical" 
massmaps (iS aha fc Winiamsll2QQ4r ). They define a phys- 
ical massmap as one that is positive everywhere and 
satisfies restrictions on maximum pixel-to-pixel varia- 
tion and direction of the massmap gradient. Additional 
constraints are optionally imposed: inversion symmetry 
about the axis and a minimum radial slope. And early 
incarnations of PixeLens found those massmap solutions 
which most closely followed the light distributio n (as 
advocated in Saha fc Williams 1997 ; Abdelsalam et al.l 
1998). 

After some experimentation, we have developed a new 
measure of physicality. Rather than imposing rigid ar- 
bitrary constraints, we assign a figure of merit to each 
massmap solution based on the following physicality 
traits: 

1. Positive mass everywhere within the convex hull 

2. Low mass scatter in each radial bin, thus preferring 
azimuthal symmetry 

3. No "tunnels": penalty for individual mass pixels 
within the convex hull being "too low" relative to 
others at similar radius 

4. Overall smoothness: minimal pixel-to-pixel varia- 
tion within the convex hull 

5. Average mass in radial bins decreases outwards 
(penalty for increasing outward) 

Below, we motivate this physicality list which has been 
chosen for the purposes of this paper. Note the user has 
the ability to modify this list rather straightforwardly 
given the freely- available LensPerfect software. 

The first trait is our only rigid constraint, obviously 
required of any physical massmap. Note that our 
massmaps are only constrained within the convex hull. 
Our basis functions do necessarily yield negative mass 
outside this region, but this is not necessarily a concern. 
We note that we are able to "get the right answer" inside 
the convex hull (Figs. [9l [TQ|) . And if negative mass lies 
in a symmetric ring outside the convex hull, then it will 
have no effect on the image positions inside. 

What we would like to avoid is one or more large iso- 
lated pockets of negative mass off to one side outside the 
convex hull. These may arise as "corrections" to a solu- 
tion which is not quite accurate within the convex hull. 
Large positive mass clumps outside the convex hull are 
similarly undesirable. Our second physicality trait serves 
to beat down these clumps, along with other positive 
mass clumps and underdense pockets within the convex 
hull. We may worry that this biases against the presence 
of subhalos. But turning this around, if we are biasing 
against subhalos, then we can be more confident of any 
subhalos that do arise in our solutions. And Occam's 
razor would dictate that the massmap solution with the 
fewest subhalos is most likely. Also note we are talking 
about subhalos in large sub-clumps. Smaller satellites 
generally cannot be resolved. (Our resolution is limited 
by the density of multiple images observed in the lens 
plane.) In the future, we may wish to perform tests to 
determine exactly how well LensPerfect recovers different 



amounts of substructure in large subhalos, given different 
numbers of multiple images. 

Underdense pockets incur extra penalty via our third 
physicality trait. While a small overdense area may be 
due to substructure, a small underdense area is much less 
physical. As our massmap is integrated a long the line of 
sight, an underdense pocket suggests that a "tunnel" has 
been carved through the entire mass distribution. This 
is not very plausible, and thus we penalize such tunnels 
more severely than over densities. 

We further beat lumps and pockets out of our solutions 
by minimizing pixel-to-pixel variation of mass (our fourth 
physicality trait). 

Finally, we penalize (without forbidding) mass profiles 
which increase with radius. Of course this might bias 
against spectacular fin dings like thering-like structure 
in CL0024 reported bv lJee et all (|2QQ7[ ). But we believe 
that it is perfectly healthy to bias against such spectacu- 
lar solutions. If a more pedestrian solution may be found 
which exhibits no ring structure, then it is likely that 
no ring structure exists. And if such a ring-like struc- 
ture were to persist in our solution despite this bias, we 
could be that much more confident of its existence. Ad- 
ditional tests could then be performed to show how well 
LensPerfect recovers ring structures in known massmaps 
and whether it might "recover" rings when they are not 
present. Finally, we note that in CL0024, the ring-like 
structure is detected (and predicted in simulations of a 
cluster collision) to lie well outside the strong lensing re- 
gion. 

Note that our physicality traits 2, 3, and 5 assume 
at least a rough azimuthal symmetry. We have yet to 
experiment with mass distributions which are strongly 
multi-modal. While our current method will suffice for 
many clusters, including Abell 1689, we will probably 
need to revise our penalty functions to analyze a highly 
asymmetric cluster such as Abell 2218. 

Even having settled on a set of physicality traits, it is 
unclear how exactly to calculate penalty functions based 
on them. And once these penalty functions are estab- 
lished, how should they be weighted relative to each 
other? The goal is for equally "offensive" massmaps to 
be penalized equally. Again this is highly subjective, 
but after much experimentation, we have devised a total 
penalty function which works well. The details are given 
in the Appendix [Al 

2.5. Massmap Optimization 

With a good penalty function in hand, we may then 
proceed to find those source positions and redshifts which 
produce the most physical massmap. Given many lensed 
galaxies, this is far from trivial. For example, 19 source 
positions (x, y) yield 38 parameters which must be op- 
timized over, plus up to 19 redshifts with their associ- 
ated uncertainties. And for each iteration, we must ob- 
tain a low resolution massmap to calculate our penalty 
function. This takes from less than a second to a few 
seconds as more images are added. Thus we must 
choose a scheme which efficiently navigates our 38(+19)- 
dimensional parameter space, minimizing the number of 
iterations necessary. Fortunately, we may use a few tricks 
to converge to good solutions fairly quickly. 

The first trick is commonplace in strong lensing 
massmap reconstruction methods: we build the massmap 



LensPerfect 



7 




Fig. 4. — Upper left: Massmap obtained using a single lensed 
galaxy (black): object #1 in our mock A1689 model. Based on 
this solution, a galaxy #2 (red) is delensed back to the source 
plane. As we see in the enlarged plot (top right) the de-lensed po- 
sitions (now plotted as circles) do not exactly align. We take the 
(weighted) average of these source positions, labeled with the (*) 
and "avg", and assign this to the source position for this galaxy. 
Also plotted (triangles) are the true source positions for objects #1 
and #2. Bottom left: Massmap solution obtained once both source 
positions have been assigned as above. By perturbing these source 
positions, we obtain a range of solutions all of which perfectly re- 
produce the observed image positions. Bottom right: The solution 
obtained when both sources are assigned to their true model po- 
sitions. Strictly speaking, this is a mock model so the scale is 
arbitrary. But if this were actually A1689, the scale would be 0'/40 
I pixel. 

one galaxy at a time. (That is, we add the multiple im- 
ages of each galaxy in turn as constraints to our model.) 
But we take this a step further, tearing our solution down 
and rebuilding it for every iteration. These rebuilds are 
extremely quick, as at each step, the deflection field solu- 
tion need only be calculated at the new image positions. 

As we add each galaxy to our model, we place its source 
position using the average of those predicted from the 
current solution (Fig. [4]). Each image will delens back 
to a slightly different position (as these images have not 
yet been constrained). By taking the average of these 
(weighted, as discussed below), we can obtain a good 
initial guess as to its source position. 

Given this source position we will obtain a new 
massmap solution. But we can also perturb this source 
position and obtain a different massmap solution instead. 
Thus we iterate and search for that source position which 
yields the most physical solution according to our penalty 
function. We perform this 2-dimensional optimization 
usmg theH owelll (Il964l ) routine^^ included in the SciPy 
Python package.^ If the redshift of this galaxy is un- 
certain, we may also optimize the redshift at this time 
(including a penalty if it drifts too far from its expected 
value). (For this one-dimensional optimization, we use 
the simple golden section search method, also included in 
SciPy.) Once the source position and redshift have been 
optimized, we may proceed to adding the next galaxy. 

Before proceeding to the next galaxy, we may choose 
to re-optimize all previous source positions and redshifts. 
This "reshuffling" can improve the overall solution. But 

More elaborate routines are available, but these may actually 
be less efficient. Specifically, calculating gradients of the penalty in 
source position space would require two extra and time-consuming 
function evaluations for every iteration. 

http://www.scipy.org 



once, say, 10 galaxies have been placed, we may worry 
that our solution has been "locked in". Attempting to 
re-optimize the third galaxy position will not be very 
effective, as it is now tightly constrained by the positions 
of the other nine galaxies. 

Thus we use a ffexible parameterization in which the 
perturbations to the predicted source positions (as de- 
scribed above) constitute our free parameters to be op- 
timized. That is, we maintain a list of source position 
offsets (one for each galaxy, initially set to (0, 0)), and it 
is this list which we optimize, rather than the source posi- 
tions themselves. Every time we rebuild our massmap so- 
lution, we obtain each source position as shown in Fig.Hl 
but then we offset it according to our list of offsets be- 
fore placing it and proceeding to add the next galaxy. 
Note that each offset affects not only the source position 
of that galaxy, but also (as this modifies the solution) 
those of all galaxies that follow. In essence, source posi- 
tions move and adjust with each other. 

As mentioned above, we determine each new source 
position by taking weighted averages of those predicted. 
We assign more weight to a delensed position for which 
the image is close to an already established image. For 
example, in Fig. [H a red image near any black image po- 
sition will have its source position given more weight in 
the average. The idea is that if the deflection field has al- 
ready been established at a (black) image position, then 
the deflection fleld at the nearby (red) point should be 
similar. We have found this weighting scheme yields bet- 
ter source positions (closer to the true model positions) 
than straight averaging. We also had some success giving 
more weight to sources with images at large radii (where 
there is less mass available to support rapid changes in 
the deflection fleld). In the end, the exact scheme is 
somewhat irrelevant as offsets from these average posi- 
tions will be perturbed and optimized. But it is possible 
that a better scheme would yield quicker convergence. 

Above we have described all the details of our optimiza- 
tion scheme, except for how we begin! That is, where do 
we put our flrst source position? Before any galaxies are 
placed, we have no massmap solution, and thus no pre- 
dictions for the flrst source position. Fortunately just 
about any choice will do for the flrst source position, 
as all yield identical (or nearly so) solutions within the 
convex hull. Even after additional sources are added, 
shifting the flrst source position has little effect on the 
overall solution. The entire source plane basically shifts 
along with it, yielding a nearly identical solution within 
the convex hull. This is a well-known lensing degener- 
acy, but certain shifts in the source plane do yield more 
physical solutions than others, as we demonstrate in § 
13.51 Correct source positions yield a solution which is 
more symmetric outside the convex hull. 

2.6. Constraining Image Fluxes and/or Shears 

Image fluxes may provide additional constraints to 
the mass model. To constrain the fluxes (and shears) 
of lensed images, previous authors have added terms 
to their matrix equation relati ng; source pos i tions, 
lens para meters, and o bservables ("Evans fc WittI l2QQ3l : 
Congdon fc KeetonI [2005: Diego et al. 2007). However 
we prefer not to interfere with our Eq. [Ml which, given 
proper basis functions, is guaranteed to obtain a per- 
fect interpolated solution of any input deflection vectors. 
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Fig. 5. — Method used to constrain image fluxes. A small box is 
constructed around each image (bottom-right), and it is delensed 
back to the source plane (top-left) using the initial massmap so- 
lution (that derived using the image positions only). The inset 
zooms in on the delensed positions (circles). The ratio of the area 
of this parallelogram to the ratio of the "box" in the source plane 
yields the magniflcation of this image. We then adjust this ratio by 
adjusting the size of the delensed image (modifying our deflection 
fleld at the corners of the box). In this case we require a larger 
magniflcation, so we decrease the size of the source image to that 
deflned by the triangles. Our image box is constrained to delens 
to this smaller source, thus adjusting its magniflcation. In practice 
we flnd that rather than adding four constraints to the deflection 
fleld, we need only constrain two points, one above or below and 
the other to the left or right of the image. 



Fig. 6. — Input massmap to be recovered with LensPerfect given 
mock lensing data. This is one of the ea rly massmap solutions of 
Abell 1689 obtained by Bro adhurst et al.l ([2005). Mass is plotted in 
units of critical density for a source at redshift inflnity. The black 
line reappears in Fig. [7] where its meaning is readily apparent. 

we need only constrain two extra points, one above or 
below and the other to the left or right of the image. 
This not only helps reduce computing time, but it also 
keeps our number of free parameters more in line with the 
number of observable constraints. We could keep the two 
numbers exactly equal by adding a single constraint. But 
this fails to properly constrain the flux, instead squeezing 
mass out to the sides as when one sits on a balloon. 



Fortunately, relative image fluxes (magnifications) and 
shears are determined by local rates of change in the 
deflection field and thus may be easily and precisely con- 
strained by adding deflection field constraints. 

In this subsection we describe how to add flux and/or 
shear constraints for multiple images. (Also see § 13 61 ) 
We do not yet have a method for incorporating weak 
lensing shear or any other constraints from singly-imaged 
galaxies (though see § 14. 3p . 

To constrain the fluxes of multiply-imaged galaxies, we 
begin by obtaining an initial massmap that reproduces 
the image positions. Next, we construct a small box 
around each multiple image and delens each back to the 
source plane. We measure the area of each delensed box 
(now a parallelogram in the source plane) and compare 
to its original area in the image plane. The ratio of these 
areas yields the magnification for that image, at least as 
predicted for our initial model. Now, by modifying the 
deflection field at the corners of the box, we can adjust 
the relative sizes of the lensed and unlensed boxes so 
that they produce the proper (observed) magnifications 
(Fig. [5]). Given these new deflection field inputs, which 
encode both position and flux constraints for the multiple 
images, we obtain a new massmap solution that perfectly 
reproduces the observed positions and fluxes. 

We note that image shears could be constrained in 
a similar manner, if desired, by adjusting the relative 
shapes of the lensed and delensed boxes. And observa- 
tional uncertainties may be incorporated by adding mul- 
tiple realizations of noise to the measuren ients and find- 
ing a perfe ct solution to each, as done bv lEvans fc WittI 
([2003 ^ and Congdon & Keeton (^2005). 

In practice we find that rather than adding four con- 
straints to the deflection field for each flux measurement. 

Of course we can also calculate the magniflcation using Eg. [28] 



2.7. Computational Efficiency 

Other strong lensing analysis methods (with the excep- 
tion of PixeLens) face a dilemma over whether to min- 
imize scatter in the source or image plane. The former 
choice, while less robust and subject to possible biases, 
is often chosen for computational efficiency. LensPerfect 
does not have to make this choice as both source and 
image positions are always perfectly constrained. 

And LensPerfect is computationally efficient. Once 
all source positions and redshifts have been established 
(along with image positions), our massmap solution co- 
efficients are obtained "instantly" (in a fraction of a sec- 
ond) via direct matrix inversion without need for itera- 
tions. Evaluating this massmap solution on a grid, while 
not quite "instant", is still very fast. On a Mac Power- 
book G4 laptop, given Ni = 93 multiple images, the 
massmap can be evaluated on a Np = 2500 = 50 x 50 
grid in 3 seconds, scaling with NiNp. 

3. APPLICATIONS 
3.1. Massmap Recovery Test 

Here we demonstrate our technique given a mock 
galaxy cluster with simulated gravitational lensing. More 
important than the ability to perfectly reproduce all mul- 
tiple image positions is the accuracy to which we can 
recover the true massmap distribution. 

We would like to test LensPerfect for a case simi- 
lar to that observed in Abell 1689. Thus we create 
a mock galaxy cluster "Babell 1689", which is very 
similar to Abell 1689. In fact, the massmap of Ba- 
bell 1689 (Fig, p ) is a ctually a solution obtained by 
iBroadhurst et al.r(|2005l ) from their analysis of Abell 
1689. But for our purposes, Babell 1689 is just a mock 
galaxy cluster. We use it to gravitationally lens 19 mock 
galaxies at redshifts between 1 and 5.5, producing 93 
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Fig. 7. — The input massmap (Fig. [6]) is used to lens 19 model 
galaxies to produce 93 multiple images at the positions shown here. 
Each image is colored and labeled according to its source galaxy. 
The source galaxies are plotted as small dots framed by a gray 
square. The black line is the "convex hull" which traces the outer- 
most images. Inside this region, the deflection fleld is interpolated 
and thus our solution is well constrained. Outside, the extrapo- 
lated deflection fleld allows a much wider range of solutions. 
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Fig. 8. — Black: Deflections of the 93 images in Fig. □ All 
deflections are normalized to a source redshift of inflnity, and the 
tail of each vector has been moved outward to the corresponding 
image position. Red: LensPerfect interpolation of the deflection 
fleld with Ro = 700 in the units plotted. The interpolation is 
curl- free and exactly matches the vectors at the given data points. 

multiple images (Fig. [7j). This is simil ar to the number 
of multiple images (106) identified by iBroadhurst et alJ 
(2005). 

We stress that we are not analyzing Abell 1689 here. 
Our analysis of Abell 1689 and its multiple images will be 
published in an upcoming paper. Here we are analyzing 
the mock cluster Babell 1689 given its mock multiple 
images. 

The 93 mock multiple image positions along with 19 
source positions and redshifts are fed into LensPerfect 
as input. Fig. [8] shows the input defiection field scaled 
to a source redshift of infinity along with a LensPerfect 
curl-free interpolation. The solution was obtained us- 
ing the Wendland function given in Equation [iTl with a 
scale factor of Ro = 700 in the units plotted. One half 
the divergence of this defiection field gives the LensPer- 
fect massmap solution (Fig. [9|). Note that it is basi- 
cally a low-resolution interpolation of the input massmap 
(Fig. [6]). No assumptions were required about the form 
of the massmap, and yet the general form and significant 
features of the input massmap are recovered. Were more 




Fig. 9. — LensPerfect massmap solution given the 93 images in 
Fig. Ul with known source positions and redshifts and setting Ro = 
700 in the units plotted. The massmap is simply half the divergence 
of the deflection fleld plotted in Fig. [8] The scale is identical to 
Fig. El the input massmap. What we obtain with LensPerfect given 
93 images is basically a low-resolution interpolation of the input 
massmap. Note that the signiflcant features are reproduced inside 
the convex hull (the black line). Outside this line, the solution is 
very poorly constrained, and can be varied simply by changing Rq. 
Here, the mass outside goes a bit negative (we have neatly clipped 
this from our colormap). But when we increase Ro to 1500, we 
flnd the mass is positive within the entire region plotted, while the 
massmap inside is barely affected. 




Fig. 10. — LensPerfect massmap solution given the 93 images 
in Fig. O with known redshifts but without knowledge of the in- 
put source positions. We found optimal source positions, or those 
which produced the most physical massmap, as judged by the cri- 
teria we describe in the text. 

than 93 input images given, the resolution would increase 
and finer details would be resolved 13.31) . 

Note that the IBroadhurst et all (|2QQ5l ) solution may 
appear to resolve very fine detail in the Abell 1689 
massmap absent from our reconstruction of Babell 1689 
presented here. But this is just a result of their assump- 
tion that some component of the Dark Matter traces 
light. LensPerfect makes no such assumption. 

In practice, we will not have knowledge of the source 
positions. So we repeat our analysis, but without pro- 
viding the source positions as input. Instead we use the 
source position optimization method outlined in §§ 12.41 
and 12.51 and detailed in the Appendix [Aj The "most 
physical" massmap solution we find is shown in Fig. [TOl 
It is very similar to that obtained with the true source 
positions (Fig. [9]). In this case the only information we 
input to LensPerfect were the image positions and red- 
shifts, which we assumed had been measured with perfect 
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Fig. 11. — Radial profiles of model (Fig. [6]) and recovered 
massmaps (Figs. [9] and I10|) . Points within the convex hull are 
binned in radius, and the mean are plotted as symbols with lines 
demarking the minima and maxima. 

14 

precision. 

Radial profiles of the three massmaps (model, recov- 
ered with true source positions, and recovered with op- 
timized source positions) are compared in Fig. [Til Only 
points within the convex hull are considered. Agreement 
of the mean mass is excellent except for some deviation 
in the center and slight deviations toward the outside of 
the convex hull. The fine structure of the mass peaks 
is not recovered in these modest-resolution (93 multiple 
image) massmaps. 

3.2. Source Position Recovery 

How well do our optimized source positions match the 
true source positions? Inaccuracies in the source posi- 
tions propagate to the deflection field and thus to our 
massmap. Modest shifts of the entire source plane are 
acceptable, however. This well-known degeneracy does 
not affect the solution inside the convex hull (see § 13. 5p . 

With this in mind, we plot (Fig. [12]) our optimized 
and true source positions from the previous subsection. 
A constant shift of a few pixels has been added to the 
optimized positions to bring them in line (on average) 
with the true positions. This is justified (not just for our 
method, but for any method) as shifts in the source plane 
are a degeneracy in the problem. A massmap solution 
obtained with one shift is virtually identical and no less 
accurate inside the convex hull than a solution obtained 
with a different shift 13. 5p . 

Applying this shift, we find our recovered source po- 
sitions are offset from the true source positions at the 
rate of 1.05 pixels on average for the 19 systems. If this 
were A1689, 1.05 pixels would translate to 0'.'42. But 
this value cannot be directly compared to any published 
results for A1689 (or any other cluster). Source position 
offsets are likely inherent to all methods, but they can 
never be measured in practice since the true source posi- 
tions are unknown. The only measurable quantity is the 
scatter of (delensed) source positions within each multi- 
ple image system. For LensPerfect, this scatter is zero. 

■"^^ Redshift uncertainties vary greatly from one study to the next. 
Thus rather than attempting to present a test which includes "typi- 
cal" redshift uncertainties, we propose that such tests be performed 
on a case- by-case basis. We note that we have had success in mod- 
eling real-life data such as the actual Abell 1689 multiple images 
complete with their redshift uncertainties (Coe et al., in prep.). 
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Fig. 12. — Optimized source positions compared to the true 
input source positions. All of the optimized positions have been 
offset by (+2.1, -6.2) to bring them in line (on average) with the 
true positions. (The solution is very insensitive to such shifts, 
especially within the convex hull.) The optimized source positions 
are slightly biased toward a solution with higher magnification. 
The solution they yield (Fig. I10|) is perfectly valid. And it is a 
simple exercise for the user to spread the source positions further 
apart and explore solutions with lower magnifications. 

In other methods, the errors due to scatter and offsets 
may be cumulative. 

Our optimal source positions are also biased toward 
a solution with slightly greater magnification than the 
true solution. This is most likely a consequence of our 
physicality measure which rewards smoother massmaps. 
A smoother massmap will have a shallower profile and 
thus higher magnification. While we have made every 
effort not to bias profile slope, some small bias may still 
remain. We will work to reduce both the offsets and this 
bias in the future. In the meantime, it is a simple exercise 
for the user to spread out the optimal source positions 
and explore solutions with lower magnifications. In fact, 
this should be a part of any comprehensive analysis. 

We must keep in mind the ultimate goal in our analysis. 
The goal is not to obtain the single best massmap, but 
rather to determine the range of massmaps which pro- 
duce reasonable solutions. Both of the massmaps pre- 
sented in § 13.11 are perfectly valid solutions which per- 
fectly reproduce the 93 multiple image positions. We 
cannot know which is more accurate without knowledge 
of the source positions. And of course this knowledge 
is unattainable in practice. This is just an inherent de- 
generacy in the problem. We are developing methods to 
thoroughly explore the solution space and will report on 
them in future work. 

3.3. 1,000 Multiple Images 

The resolution of a LensPerfect massmap is dictated by 
the density of multiple images detected. Each multiple 
image samples the deflection field at a given location. 
In the gaps between these images, the deflection field 
must be interpolated, and the exact form of the massmap 
becomes less certain. 

To demonstrate this, we consider the massmap solution 
we may obtain given 1,000 multiple images. Rather than 
performing mock lensing to produce multiple images as 
in the previous subsection, here we will simply sample the 
deflection field at 1,000 points. We restrict these samples 
to an area of similar size as before (a circle of radius 150 
pixels). And we ensure that all samples are separated 
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strate the "steepness" or "mass-sheet" degeneracy 
(o riginany dubbed the " magnification transformation" 
in iGorenstein et alJ Il988l ) , which states that the mass 
everywhere may be replaced by 



Fig. 13. — Massmap solution to that shown in Fig. [6]given the 
deflection fleld sampled at 1,000 points within the area outlined in 
black. Fine detail in the input massmap is faithfully reproduced. 
While this demonstration is idealized, we may expect a massmap of 
similar quality given a lens which produces 1,000 multiple images. 



by 2 pixels or more. Given these 1,000 samples, we then 
interpolate to find the defiection field elsewhere. Our 
massmap solution is shown in Fig. [131 Comparing with 
the input massmap (Fig. [6|), we can see that very fine 
detail is faithfully resolved. 

This setup is a bit idealized. The "multiple images" are 
spread fairly uniformly (randomly) about a circle within 
the field. In practice we are more likely to find clumps 
and voids of multiple images due to the magnification 
pattern, obscuring foreground galaxies, and physically 
linked background galaxies. Also, we have assumed the 
equivalent to perfect knowledge of the source positions. 
Nevertheless, this massmap gives us a rough idea of the 
resolution we may expect given observations deep enough 
to reveal 1,000 multiple images. 

In practice, we would need to find and optimize a large 
number of source positions. 1,000 multiple images may 
be produced by about 300 background galaxies. This 
saddles us with 600 free parameters. But even in our 
93-image system, we find that source positions added to- 
ward the end are already fairly well constrained by those 
added previously. So we can speculate that the next 
200+ source positions may similarly "fall into place", 
making the computational challenge more manageable. 

3.4. Fewer Constraints and the Role of Rq 

While it is useful that LensPerfect can produce such 
detailed massmaps when given so many multiple images, 
what happens when LensPerfect is given far fewer con- 
straints, say a single quadruply-imaged galaxy? Does the 
sum of four of the oddly-shaped basis functions depicted 
in Fig. [3] yield a reasonable massmap solution? 

Yes it does. Fig. [14] (left-hand panel) shows the 
massmap solution obtained given a simple symmetric 
"Einstein cross" 4-image configuration. The solution is 
exactly the same (to within Hi = 0.01) if we form a more 
complete "Einstein ring" with 16 image constraints as 
shown. But this is just one possible solution, obtained 
with Ro equal to ten times the Einstein Radius. We 
begin to explore other solutions by varying the width Ro 
of our basis function (Eq. [T7|) . Radial profiles of these 
different solutions are shown in the right-hand panel. 
By changing this single parameter, we neatly demon- 



1- K.' = X{1- K.), 
a^' = Aa^ + (1 - A) 



(31) 
(32) 



without affecting the observed image positions (unless 
multiple galaxies of different redshifts have been lensed). 
Note that the new mass hz^ is steeper than the previous 
mass hzhy di factor of A. 

The different LensPerfect solutions follow the transfor- 
mation in Eq. [31] very closely, although not exactly. We 
have rescaled the green Ro = 100 curve via this transfor- 
mation such that its peak aligns with the Ro = 55 pro- 
file. It now aligns with the full i^o = 55 curve very well 
within the Einstein radius, but not perfectly. (Note that 
the Einstein radius is also the convex hull in this case. 
Thus we do not expect to be able to constrain the solu- 
tion outside.) Also note in the zoomed subplot, that the 
different profiles attain k = 1 dit nearly the same radius 
R, but not exactly, as they would if they had resulted 
from the transformation in Eq. [311 Thus, the degener- 
acy we probe by simply varying Ro is very similar to but 
slightly different than the classic and simplest form of 
the "mass-sheet degeneracy" . 

We note that these solutions all have fiat i<i{R) slopes 
at the center, corresponding to a fiat slope in p(r), the 
three-dimensional density, as well. This slope gradually 
increases outward from the center. This i s the same be- 
havio r seen in many CDM sim ulations (Navarro e t al.l 
120041 : [Merritt eTall l2005l. [20061. an d references therein, 
although see iDiemand et al. IM), which suggest that 
Dark Matter halos may attain the same "Persic (^9681 ) 
density profiles that we observe in the light profiles of 
elliptical galaxies. Thus the LensPerfect mass profiles 
obtained when given only four image constraints do ap- 
pear to be reasonable. But we stress that these solutions 
are far from unique, and many other mass profiles are 
possible. 

3.5. Source Plane Shifts 

Theory tells us that given a single source or even mul- 
tiple sources at the same redshift, the entire source plane 
can be shifted, and a new massmap solution can be found 
which does not affect the image positions. The degener- 
ac y was originally named t he "prismatic transformation" 
bv IGorenstein et all (|l988l ). who explained that this shift 
can be accomplished by adding an infinitely long and 
thin mass stick to the solution off to one side of the im- 
ages. Of course such mass sticks are not very physical. 
Thus, we can overcome this degeneracy by simply find- 
ing that solution which does not resort to the addition 
of mass sticks! This is why we "know" that Einstein 
rings are produced by an on-axis alignment of the lens 
and source galaxies. The on-axis configuration requires 
a simple symmetric massmap, while any off- axis source 



As noted by S aha & Williams' (2006), the term "mass-sheet" 
degeneracy may lead to the mistaken belief that a constant mass 
sheet is what may be added without affecting the images. In fact 
the mass must also be rescaled, or multiplied, in concert. 
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Fig. 14. — Left: LensPerfect massmap solution with Rq = 100 given an Einstein cross: 4 multiple images located at the black points, 
and the source position at the center. Given a more complete Einstein Ring (images at all 16 black and grey points), the massmap solution 
is identical to within k = 0.01. The circle plotted inside the Einstein radius is the = 1 contour. By changing Rq (the width of our basis 
function), we obtain solutions with different radial profiles, the range of which is demonstrated by four examples plotted at right. This 
demonstrates the classic "steepness" or "mass-sheet" degeneracy 1 — hi^ = X(l — almost exactly. The dashed green line is the Ro = 100 
solution rescaled via this transformation to match the peak value of the Ro = 55 solution. The two line up nearly perfectly, indicating that 
the degeneracies shown here are very similar to but not exactly the classic form. Also note in the zoom inset, that the different solutions 
have = 1 at slightly different radii, another indication that the classic mass-sheet degeneracy relation is being followed approximately, 
but not exactly. 
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Fig. 15. — Einstein Ring produced by a galaxy offset from the 
ring. While this configuration is possible, it is highly unlikely that 
an offset source galaxy and asymmetric lensing mass distribution 
would conspire so precisely to produce an Einstein ring in Nature. 
But do note that the massmap solution inside the ring is nearly 
identical to that obtained with the source position placed at the 
center of the ring. 

position would require a less physical massmap, perhaps 
requiring the addition of a mass stick. 

We say "perhaps" because LensPerfect is able to find 
solutions to such off-axis Einstein rings without resorting 
to mass sticks. (Fortunately LensPerfect has no mass- 
stick basis functions to work with.) One such solution 
is shown in Fig. [151 This mass distribution, while more 
plausible than a mass stick, is nevertheless highly im- 
probable. The proper amount of mass must be added 
off-axis in the correct configuration to refocus the light 
just enough toward the center of the lens where it then 
produces the Einstein Ring. It is very unlikely that such 
a precisely-tuned mass configuration would occur in Na- 
ture. Much more likely is the on- axis alignment, in which 
case any simple symmetric mass configuration will do. 
Our optimization procedure (§§ \2.4\ I2.5p quickly con- 
verges to the on-axis source position as being most likely, 
as off-axis solutions are penalized for being asymmetric. 

But we should not be surprised to find modest source 
position shifts in our solution either (§ 13. 2p . The mass 
within the convex hull is extremely insensitive to these 



shifts. (A constant shift in the source plane does not 
affect the divergence of the defiection field at the multiple 
image positions. Thus the massmap within the enclosed 
region experiences only negligible changes.) 

3.6. Constraining Image Shapes and Sizes 

Additional constraints may be derived from gravita- 
tionally lensed images by requiring the mass model to 
correctly reproduce not only image positions, but also 
their sizes, shapes, orientations, and fiuxes. A method 
to constrain the fiuxes (and shears) of unresolved images 
was discussed in § 12.61 But observed fiuxes (and espe- 
cially shears) may be uncertain. More precise constraints 
may be obtained when the multiple images are resolved 
and extra knots may be identified within them. 

Fig. [16] shows two of the four multiple images of a 
z = 4.92 galaxy pr oduced by the galaxy cluster MS 1358 
(|Franx et al.l ll997l ). Six knots are identified as common 
in both of these images (B & C), with three of the knots 
also being identified in image A (not shown). For each 
knot in turn, all of its images are constrained to originate 
from the same source position. The defiection model is 
updated after each knot is added. The final defiection 
solution yields the source plane images in the right-hand 
panel (see also Fig. [T7j). Note the identical alignment of 
the two images (B & C). At this point, the de- lensed im- 
ages may be co-added t o obtain greater depth, if desired 
re.g.. "Colle v d~ZI [T996). 

But this capability allows us to do more than simply 
produce pretty de-lensed images. The extra constraints 
provided by image knots can add significant detail to the 
massmap solution. Fig. [TSl compares (left) the MS1358 
massmap solution obtained when using only a single im- 
age position (the brightest knot) for each object and 
(right) the solution obtained when constraining all knot 
positions. In the latter, the mass is more stretched to- 
ward a second cluster galaxy which proves crucial to tra- 
ditional parametric mass modeling. 

4. DISCUSSION 
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Fig. 16. — Left: Two images of a quadruply-lensed z = 4.92 galaxy behind MS1358. Six matching knots are identified and labeled in 
each of these two images. A cluster galaxy has been subtracted from this BVi' z' JH color image. Right: The same image de-lensed ba ck t o 
the source plane given the LensPerfect massmap solution. All knots have been constrained to align in the source image. (Also see Fig. I17h 




Fig. 17. — Images B and C de-lensed back to the source plane using LensPerfect. The six knots have been constrained to align exactly 
in the de-lensed images. 




Fig. 18. — MS1358 massmap solutions obtained by LensPerfect. Left: Only the positions of the brightest knots are constrained. Right: 
The positions of six matching knots are constrained in images B and C along with three knots in image A and one in image D. Note that 
the colormaps in the two plots do not have the same scale. 
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We now describe the current limitations of the 
LensPerfect method and plans for future capabilities 
LensPerfect may provide us with. We also discuss the rel- 
ative merits of both "parametric" and "non-parametric" 
methods, elaborating on points made in the introduction. 

4.1. Imperfect Perfect Solutions 

The "perfect" in LensPerfect refers to the reproduc- 
tion of multiple image positions attained by its massmap 
solutions. But these solutions may be perfectly wrong 
and even unphysical if incorrect source positions and/or 
redshifts (and/or multiple image identifications) are pro- 
vided. We have developed a method to find source posi- 
tions and redshifts which produce reasonable solutions. 
We do not claim this method is perfect; it appears to 
work well but may be improved upon in the future. 

LensPerfect massmap solutions do not guarantee 
against predictions of additional multiple images which 
don't exist. This is an issue common to many methods, 
and it has been argued alternately that this should ei- 
ther more or less be a concern. The argument against 
its importance is that the prediction of additional multi- 
ple images is very sensitive to local substructure. Slight 
modifications might be made to this substructure, it is 
argued, to make these extra predicted multiple images 
disappear. LensPerfect may finally give us a tool capa- 
ble of easily demonstrating this. We may be able to add 
defiection field constraints that tweak the massmap just 
enough to eliminate unwanted images. This idea has not 
yet been tested. In the meantime, we will simply check 
each solution for the correct number of multiple images. 
We find our models generally do not produce extra mul- 
tiple images, as long as reasonable source positions are 
assumed. 

4.2. Imperfect Solutions Superior? 

Would it help to loosen the restriction that our 
massmap solutions perfectly reproduce all the multiple 
image positions? That is, might we obtain more accu- 
rate solutions by allowing our solutions to be imperfect? 
There are two parts to this question. 

First we consider measurement uncertainties. We have 
already discussed how we deal with redshift uncertain- 
ties in § 12.51 We allow the redshift to vary within 
the uncertainties as we optimize the massmap. The 
same could be done for other uncertain measurements. 
Image position uncertainties are generally very small 
(about one pixel). But when adding fiux / shear con- 
straints, we should certainly include the corresponding 
unce rtainties in our optimization. Evans & Witt (2003) 
and Con gdon fc KeetonI (|2005i) followed this approach in 
their studies of galaxy lenses. By including uncertainties, 
they noted modest improvement in their solutions. 

The second part of the question is more fundamen- 
tal. Could a less perfect massmap be more accurate? 
This does prove true in the SLAP method (Diego et al. 
[2005a). They purposely leave some scatter in the source 
plane to avoid solutions which are "biased" with "a lot 
of substructure". But this is likely a result of their for- 
mulation of the problem. They consider all pixels in each 
lensed image and find that massmap solution which min- 
imizes the sizes of all the delensed images. Of course the 
delensed sizes should not be zero, so the delensed image 



pixels are allowed some "scatter". (They make no at- 
tempt to match internal features as we demonstrate in § 
13.61 ) Meanwhile, the proponents of PixeLens do not re- 
port such problems with their solutions which perfectly 
reproduce all image positions (Saha et an i2006'). 

We performed a test to see what an imperfect solu- 
tion would look like and if it might be more accurate. 
We began with our perfect solution to the 93 multiple 
images in § 13.11 with optimized source positions. We 
then "reoptimized" the source positions (based on our 
physicality criteria), one-by-one for each individual im- 
age. Source positions were no longer constrained to a 
single point for each multiple image system. In practice, 
they drifted by an average of 0.4 pixels. Our resulting 
imperfect massmap was not extremely different from our 
perfect massmap but not any more accurate either. Some 
of the substructure features in the center were erased as 
the physicality measure guided the optimization toward 
a smoother solution. The radial profile did not signifi- 
cantly improve nor deteriorate in accuracy. 

We believe the "1,000 points of light" massmap in 
Fig.[T3]is a convincing demonstration that perfect is best. 
When the defiection field is sampled at 1,000 points, we 
are able to map fine detail in the massmap. This fine 
structure is encoded in the exact positions of the mul- 
tiple images. Allowing for an imperfect solution would 
erase this fine structure and waste information obtained 
in the observations. 

4.3. Weak Lensing, Extended Images, and Time Delays 

It may be possible to directly incorporate measure- 
ments of both "weak" lensing and not-so-weak shear into 
the LensPerfect method. Lacking multiple identifiable 
image knots, we can constrain galaxy shapes using the 
same method we use to constrain fiuxes (§ 12. 6p . Instead 
of altering the relative sizes of the source and lensed re- 
gions, we can alter the shear. Thus we could constrain 
all delensed images to be round, or perhaps round with 
some scatter of ellipticity as measured in large scale sur- 
veys. We do not develop this idea further here. Such 
a technique would require three constraints per galaxy, 
which, given 100 galaxies, would require significant com- 
puting time (although at this stage iterations may not 
be required). 

A more practical idea may be to compare measured 
shears to those predicted from our mass models and in- 
clude the disagreement as a penalty in our optimization 
procedure. But note that this would only work well for 
galaxies within our convex hull as our solution (and thus 
shears) would be very poorly constrained outside this re- 
gion. 

Beyond simple constraints of fiux, shear, and im- 
age knots, all of the pixel information in each lensed 
image may be utili z ed in the model derivation (e.g;., 
Warren fc D^l2003l : ISuvu et al.ll2006l : iKoop mans et al.l 
20061 ). This is especially desirable in strong galaxy- 
quasar lensing, in which observable constraints are less 
plentiful than in cluster lensing. It is unclear how to im- 
plement this in LensPerfect, except as another penalty 
in our optimization routine. 

Another constraint often used in strong galaxy-quasar 
lensing studies is time delays observed among different 
images. Time delays are not local functions of the defiec- 
tion field, and as such, we cannot constrain them directly 
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with LensPerfect as we constrain fluxes and shears. But 
for a given solution, they can be calculated directly from 
derivatives of our basis functions and compared to those 
observed, yielding another penalty for our optimization 
routine. 

4.4. 'Tarametric^^ vs. Non-parametric^^ Methods 

Methods are often classified as "parametric" or "non- 
parametric", depending on whether a clear physical 
parameterization is used to construct the proposed 
massmap solutions. "Grid-based" methods are often 
referred to as "non-parametric", even though strictly 
speaking they do have parameters, namely the mass at 
every pixel on a grid. The real distinction to be made 
here is between "model-based" and "model-free" meth- 
ods. The former construct mass halos as physical an- 
alytical forms, while the latter do not. To give exam- 
ples, the fully model-based strong lensing analyses of 
the A1689 ACS iniages published to date have been 



tne Aibob) ALyD im ages p ublisned to date nave been 
iHalkola et al." ("20 061. 12007 ) and L imousin et"al] (l2006l). 
Mean while. Diego et alTj 200 5b). ISaha et all ([2006.) . 
and Leonard et al.l (|2007l ) have obt ained mo del- free 
massmaps based on these images. The iBroadhurst et al.l 
|2005) and Zekse r et all (|2006r ) analyses included both 
model-based and model- free elements. 

LensPerfect, while clearly parametric, is also model- 
free. The LensPerfect solutions are given as sums of ba- 
sis functions. But these basis functions have no physical 
interpretation. And this functional form is practically in- 
discernible in the final solutions. In fact, the basis func- 
tion coefficients are generally many orders of magnitude 
greater than the amplitude of the massmap. These large 
mass components all cancel out nearly perfectly with the 
"residuals" being the massmap solution. 

But while LensPerfect is model- free, it does not have 
the large number of free parameters typical of grid-based 
methods. In fact when fitting image positions only (in- 
cluding extra knots, but not fluxes), the number of free 
parameters solved for (the coefficients) is exactly equal to 
the number of constraints. Note that the source positions 
and redshifts are not solved for directly, but rather must 
be provided as input. Each flux measurement provides 
a single constraint but two free parameters are required 
to constrain it in our models. And a shear measurement 
provides two constraints, which can be reproduced with 
an equal number of (two) free parameters. 

But which are superior overall, model-based or model- 
free methods? Individual researchers may have a decided 
preference for one over the other, but in fact both types of 
methods have their strengths and each serves a purpose. 

Physical model-based massmaps allow us to test 
whether Dark Matter is distributed in certain ways, ac- 
cording to our assumptions. In principle, they allow for a 
more straightforward determination of meaningful phys- 
ical parameters. Cluster mass models attempt to simul- 
taneously constrain the forms of both the overall cluster 
halo and the individual galaxy halos (e.g.. iHalkola et al] 
|2007[ ). But degeneracies are strong between the two ad- 
ditive components. 

On the other hand, model-free massmap reconstruction 
methods allow us to directly test for the presence of Dark 
Matter free of assumptions about its distribution. In par- 
ticular, these methods make no assumptions about mass 
following light. Some of our most important discoveries 



about Dark Matter have come and are expected to come 
from cases in which mass does not follow light. Mass 
peaks offset from light peaks can provide constraints on 
the collisional nat ure of Dark Matter pa r ticles, as in 
the Bullet Cluster (Markevitch et al."2004l: IClowe et all 
[200 6) . And observations of Dark substructure (not as- 
sociated with light) may someday vindicate CDM simu- 
lations w hich predict much more halo substructure than 
is visible (|Diemand et al.ll2007l : IStrigari et al.ll2007l . and 
references therein). 

Model-free methods are also able to explore a wider 
range of possible massmap solutions, including (with the 
advent of LensPerfect) those which perfectly reproduce 
all 100+ multiple image positions. Of course, exploring 
the full range of model-free solutions can be a compu- 
tational challenge. And it is not entirely clear how to 
sort the physical solutions from those which are "less" 
physical or aphysical. 

Meanwhile, too much model flexibility has been cited 
a potential problem, as a flexible mode l may fit incor- 
rect d ata without any alarms sounding. iLimousin et al.l 
(|2Q06() claim tha t their model-based method and that 
of IHalkola etHI ([2006h were unable to fit some mul- 
tiple image syste n is inc orrectly identified in the initial 
IBroadhurs t et all (|2005l ) work, which was a bit more 
model-free. But LensPerfect, while more flexible still, 
may actually be more unforgiving than all previous meth- 
ods when it comes to incorrect multiple image identifi- 
cations. The reason is simple. "Imperfect" massmap 
reconstructions experience some offset in each and every 
predicted image position. Thus a misidentified multi- 
ple image set may have a larger than average, but 
this may be more easily dismissed in the analysis. In 
LensPerfect, however, each multiple image puts a rigid 
constraint on the deflection field. Thus a misidentified 
multiple image is more likely to cause the deflection field 
to get tangled, leading to a "less physical" (if not aphys- 
ical) solution. We stay alert for such ill-fitting multiple 
image systems as we add each to our models. And, of 
course, we must take care to avoid misidentifying multi- 
ple images from the start whenever possible. 

Another obje ction to model flexibility was raised by 
iKochane^ (|2004 ) . He ta kes exception to the a bility of the 
model-based but flexible lEvans fc WittI (|2003[ ) method to 
produce a mass model that accounts for the flux anomaly 
observed in the strong galaxy-quasar lens Q2237+0305 
even though this flux anomaly has since been shown 
to have been due to a microlensing event, which has 
now pass ed! We would ar^ue that the mass model pro- 
posed bv ' Evans fc WittI (|200 3') is but one possible solu- 
tion among several to the observed flux anomaly. All 
possible solutions should be considered, and in this case, 
microlensing is proven to be the true solution. 

In analyses of galaxy-quasar lensing, model-based 
methods enjoy even greater appeal than in cluster lens- 
ing. Perfect solutions are much easier to come by when 
there are fewer constraints (4 image positions, for exam- 
ple). And with so few image constraints, a large range 
of solutions is possible. We demonstrated one way of 
probing this solution space in § 13. 4[ and the PixeLens 
method provides another. But many choose to slash the 
solution space by making well-motivated model-based as- 
sumptions (e.g., an isothermal proflle). 

A powerful technique would be to combine LensPerfect 
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with a model-based method. We could find good (imper- 
fect) model-based solutions which leave offsets between 
the observed and predicted image positions, and then 
perfect these solutions with LensPerfect. The defiection 
offsets may be interpolated, and this "offset solution" 
added to the imperfect solution to create a perfect solu- 
tion. We implemented this idea, but the results appeared 
unruly. A more creative approach will be required if this 
capability is to be realized. 

We mention one other advantage to our method. That 
advantage is ease of use. If a simple mass model is re- 
quired, LensPerfect 's speed is hard to beat. Given a sin- 
gle lensed galaxy, LensPerfect instantly obtains a per- 
fect solution "out of the box" . And if instead given 30+ 
lensed galaxies as in Abell 1689, LensPerfect still obtains 
solutions with minimal user input. Parametric methods 
instead require the user to develop complex models ca- 
pable of fitting so many multiple images well. Previ- 
ous studies identified and measured properties of many 
cluster galaxies in order to model a "galaxy component" 
which would be added to a separate halo component in 
their solutions. LensPerfect instead bypasses this param- 
eterization process obtaining a detailed massmap free of 
strong assumptions. 

5. SUMMARY AND FUTURE WORK 

We have presented a new approach to gravitational lens 
massmap reconstruction. Given image positions, source 
positions, and redshifts, a new mathematical technique 
is used to interpolate the defiection field via direct inver- 
sion. The resulting massmap (simply half the divergence 
of the defiection field) perfectly reproduces all of the ob- 
served image positions. 

In practice, source positions are unknown. We have 
devised a method that efficiently optimizes over different 
possible configurations of the source galaxies. Each con- 
figuration produces a different massmap solution which 
is evaluated for "physicality" based on criteria developed 
as part of this work. Our criteria make only minimal as- 
sumptions about the form of the massmap. Specifically, 
they make no assumptions about the slope of the radial 
profile nor mass following light. 

We demonstrated our method on mock gravitational 
lensing data. A known massmap was used to lens 19 



mock galaxies to produce 93 multiple images. Our 
massmap solution perfectly reproduces all the multiple 
image positions while accurately recovering the known 
lens mass distribution to a resolution limited by the den- 
sity of the multiple images. We demonstrate the im- 
proved accuracy and fine detail we may expect from a 
massmap derived from 1,000 multiple images. 

We also presented LensPerfect solutions based on far 
fewer constraints, such as four multiple images of a single 
galaxy. Image fiuxes or extra knots may provide addi- 
tional constraints which can also be perfectly fit by our 
models. The number of free parameters solved for is kept 
nearly equal to the number of observable constraints. 

The LensPerfect software is easy to use and is made 
publicly available at our website (see §[1]). 

In subsequent papers, we will apply this method to 
Abell 1689 and other galaxy clusters. We will also at- 
tempt to resolve some of the "fiux anomalies" observed 
in galaxy lenses. 

In § m we discussed some possible improvements and 
extensions of our method, including the incorporation of 
weak lensing data. We also hope to better quantify the 
accuracy of our massmap recovery and compare our per- 
formance to that of other methods. The uncertanities in 
our solutions should be well determined by exploring the 
full range of physical solutions. We will also experiment 
with other radial basis functions and other recipes for 
measuring massmap physicality. 
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APPENDIX 

A. MASSMAP PHYSICALITY PENALTY FUNCTION 

To calculate our penalty function, we first evaluate the massmap solution on a 41 x 41 grid. This proves large enough 
to measure physicality, while taking less than a second to calculate given as many as 50 image positions. (As this 
evaluation is to be part of our iterative procedure to determine source positions, it should be kept as quick as possible.) 
Based on this 41 x 41 massmap, we calculate the total penalty as follows: 

1. The sum of all negative pixels inside the convex hull is multiplied by —100. (We could simply assign a penalty 
of infinity to negative massmaps, but this would create a large region of constant penalty in source coordinate 
space. Our finite and varying penalty serves better during the optimization to corral the source positions back 
toward those which yield positive solutions.) 

2. The mean mass {tz) and RMS scatter Ahz are measured in radial bins of 80 points each. The RMS scatter is 
totaled and divided by 10. Within the convex hull, the mean and scatter are recalculated in radial bins of 40 
points each. This inner scatter is totaled and multiplied by 4. Outside the convex hull, we rebin the massmap, 
yet again, in bins of 80 points each. Rather than the RMS mass scatter, this time we penalize the peak-to-peak 
variation (that is, the maximum minus the minimum) in each bin. These variations are totaled and multiplied 
by 5 and divided by the number of bins. 
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3. Any outward increase in {k) over R is measured and multiplied by 10. This is calculated both for all points and 
then again for those within the convex hull. Both penalties are added. 

4. Given the mean mass (hi) and RMS scatter Ak, calculated within the convex hull (see 2 above), we interpolate 
the l-cr lower limit {k,) - Ak, for all points within the convex hull. Deviation below this limit is measured and 
divided by 2. Note that we certainly expect some points to fall below the 1-a lower bound. (By definition, 16% 
of the points should fall below.) But the idea is to minimize these deviations. A large deviation below indicates 
a "tunnel" or unphysical dip in the surface density. 

5. Finally, we put a premium on smooth solutions within the convex hull as being the most likely. Numerical 
derivatives of the massmap are calculated at each pixel (x, y): dn/dx = \k{x + 1,^) — k,{x — 1,^)|, dhc/dy = 
|/^(x,^ + l) — K{x,y — 1)1 The absolute values of these are totaled the entire sum and divided by 5. 

The above penalties and their respective weights were defined after much trial and error. We find that massmaps 
with lowest total penalty defined as above appear to be the best behaved, or most physical. Again, other weights and 
penalty schemes are certainly possible. 
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